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We study the nonequilibrium dynamics of a many-body bosonic system on a lattice, subject to 
driving and dissipation. The time-evolution is described by a master equation, which we treat within 
a generalized Gutzwiller mean field approximation for density matrices. The dissipative processes 
are engineered such that the system, in the absence of interaction between the bosons, is driven into 
a homogeneous steady state with off-diagonal long range order. We investigate how the coherent 
interaction affects qualitatively the properties of the steady state of the system and derive a nonequi- 
librium phase diagram featuring a phase transition into a steady state without long range order. The 
phase diagram exhibits also an extended domain where an instability of the homogeneous steady 
state gives rise to a persistent density pattern with spontaneously broken translational symmetry. In 
the limit of small particle density, we provide a precise analytical description of the time-evolution 
during the instability. Moreover, we investigate the transient following a quantum quench of the 
dissipative processes and we elucidate the prominent role played by collective topological variables 
in this regime. 
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I. INTRODUCTION 

One major challenge in nonequilibrium many-body 
physics is to identify situations with a sufficient degree 
of universality, i.e. phenomena that occur independently 
of a precise microscopic realization and are largely insen- 
sitive to the specific choice of initial conditions. Several 
situations have been discussed in the literature, where 
many-body systems evolve in time towards nonequilib- 
rium steady states. The most prominent example in con- 
densed matter is certainly the electron gas exposed to a 
bias voltage pQ. In this context, truly many-body prop- 
erties, such as the effect of nonequilibrium conditions on 
quantum critical points, have been investigated [2 j. Fur- 
ther implementations of nonequilibrium many-body sys- 
tems are discussed in the context of exciton-polariton 
Bose-Einstein condensates [3 , or more recently using 
driven noisy systems of trapped ions or dipolar atomic 
gases [4 . As far as the time evolution of many-body 
systems is concerned, for closed systems we have seen a 
plethora of studies of quench dynamics [SJ [6] , thermaliza- 
tion [7, 8], and pre-thermalization [9 , as well as trans- 
port [10] . Further dynamical studies involve situations of 
crossing quantum critical points in a finite time, in the 
spirit of the Kibble- Zurek mechanism [TTJ [12], or investi- 
gations of the many-body Landau-Zener effect [T3] . 

The quantum optics toolbox for the manipulation of 
cold atomic systems [14 not only offers the possibility 
to tailor the coherent dynamics of a closed many-body 
system by Hamiltonian engineering, but also to imple- 
ment the dissipative dynamics generated by a Liouville 
operator [15 . Of particular interest are Liouville op- 
erators characterized by the existence of a unique dark 
state, i.e. a sink state whereto the system is asymptot- 
ically driven, independently of the initial configuration. 
In this case, it has been shown theoretically that the 
action of the sole dissipation allows to reliably prepare 



the system in peculiar states with interesting quantum 
mechanical features, such as phase coherence in bosonic 
systems [15], entanglement in spin systems [16] [17], or 
delocalized pairing with arbitrary symmetry in fermionic 
systems [18]. Since a dark state can be reached by dis- 
sipative dynamics but not left, such processes featuring 
a dark states clearly violate the principle of detailed bal- 
ance and hence may drive the system into a steady state 
quite different from thermodynamic equilibrium. The ex- 
istence of a dark state with known properties nevertheless 
promises a sufficiently universal nonequilibrium situation 
with a well-defined steady state. 

Furthermore, since dissipation needs not be a small 
perturbation, but on the contrary can be made the dom- 
inant contribution to the dynamics, it is natural to in- 
vestigate the effects of the interplay between unitary and 
dissipative dynamics. In thermodynamic equilibrium, 
the competition of two non-commuting microscopic op- 
erators leads to a phase transition if the ground states 
have different symmetries when one or the other oper- 
ator dominates [19]. A seminal example in the context 
of cold atoms in optical lattices is the superfluid-Mott 
insulator transition in the Bose-Hubbard model [2D] . By 
analogy with equilibrium physics, one may similarly ex- 
pect a phase transition in the steady state of the system, 
resulting from the competition of unitary and dissipative 
dynamics. 

In this paper, we demonstrate a phase transition in 
the steady state of a driven-dissipative bosonic system 
on a lattice and discuss its properties in detail. We point 
out several features of the phase transition that make it 
a genuine nonequilibrium phenomenon, not fitting into 
the well-developed theory of equilibrium quantum phase 
transitions. However, we recognize some properties that 
may hold universally, i.e. that may be shared by entire 
classes of similar systems, in the nonequilibrium steady 
state or in the time evolution. The profound difference 
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between the phenomena that we consider here and the 
equilibrium quantum phase transition manifests itself, at 
a purely technical level, in the absence of a known vari- 
ational principle for the steady state. The lack of such 
formal tool makes it necessary to evaluate the real-time 
evolution of the system, which is rather impractical in the 
many-body case. Here we develop a mean field method 
to treat bosonic lattice systems, subject to unitary and 
dissipative dynamics, which generalizes the Gutzwiller 
ansatz and allows us to describe mixed states. Within 
this approximation scheme, we deal with an effective 
master equation that acts on a reduced Hilbert space, 
is nonlinear in the quantum state, and is amenable of a 
workable numerical solution. Moreover, we identify an 
analytically tractable low-density limit where the most 
relevant information on the quantum state is encoded in 
a few correlation functions, which obey coupled equa- 
tions of motion. The solution of the equations of motion, 
which we provide for various dynamical regimes in this 
approximation, explains qualitatively the behavior of the 
system at arbitrary density. The method may also be 
used for semi-quantitative estimates in other problems 
where unitary and dissipative dynamics appear on equal 
footing. 



The paper is organized as follows. In Sec. [II] we in- 
troduce the microscopic model and provide a heuristic 
argument to justify the existence of a nonequilibrium 



phase transition in the model. In Sec. Ill we introduce 
the generalized Gutzwiller mean field approximation, dis- 
cuss the most prominent properties of the homogeneous 
phase diagram resulting from it, and point out profound 
differences with respect to equilibrium phase transitions 
in dissipative environments [21]. Sec. |IV| specifies the 
Gutzwiller approximation to the case of low particle den- 
sity and complements the numerical findings of Sec. [TTT] 
with analytical results. Sec. [V] is devoted to the analy- 
sis of a dynamical instability that arises spontaneously 
in the nonequilibrium phase diagram. This completes 
the discussion of the steady state of the driven dissipa- 



tive system. In Sec. [VI] we focus on aspects of nonlinear 
dissipative dynamics that occur in mesoscopic settings. 
We draw our conclusions in Sec. |VII| The appendices 
[AUD] provide technical details on the derivation of the 
mean field theory and on the numerical procedures used 
to solve the resulting nonlinear equations. 



This paper refines the findings that we already pre- 
sented in Ref. [22]. Here we provide a more detailed 
presentation of the material, including the discussion of 
the complete equations of motion and the explanation of 
numerical techniques. New results include the exact so- 
lution of the dynamics for a small lattice (Sec. [il]), the 
solution of the dynamically u nstab le equations of motion 
beyond linear response (Sec. VB), as well as the inves- 
tigation of the nonli near dissipative dynamics following 
phase quenches (Sec. |VI|). 



II. THE MODEL FOR COMPETING UNITARY 
AND DISSIPATIVE DYNAMICS 

A. The model 

In this work we consider a many-body system de- 
scribed by a master equation 



d tP (t) = -i[U,p(t)}+C[p(t)} . 



(1) 



The master equation defines the time evolution of the 
density matrix p(t) of an ensemble of bosonic atoms in 
an optical lattice. It consists of a unitary part, described 
by the Hamiltonian and a dissipative part, represented 
by the Liouvillian C. The Liouville operator arises from 
the coupling of the system to a bath, which is eliminated 
in the Markov approximation to yield the effective evo- 
lution 0. The validity of the Markov approximation, 
which neglects retardation effects in the bath, is rooted in 
the possibility of having strong separation of time scales 
for such quantum optical systems. We emphasize that 
this gives rise to a microscopic model which is local in 
time. Therewith, the possibility of an accurate micro- 
scopic modeling for many-body systems in terms of a 
few parameters - one of the key features of cold atomic 
many-body systems - is extended from the Hamilton op- 
erator to the dissipative dynamics as well. The unitary 
dynamics for the bosonic atoms in a <i-dimensional lattice 
is described by the familiar Bose-Hubbard Hamiltonian 
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where J is the hopping amplitude, U > the i nterac tion 
strength, and \i the chemical potential (see Sec. |IIIB and 
App. [C] for a discussion of its role in the time-evolution 
of the system). The operator bg (b\) annihilates (creates) 
a boson in the ^th lattice site and (££') indicates that the 
summation is restricted to nearest neighboring sites. {I 
is a vector of integers with the same dimension as the 
lattice.) The Liouvillian C that we consider here was 
derived in Ref. [15]. Written in Lindblad form, it reads 



C\p(t)] = \* E^«' 



(3) 



where n is the dissipation rate and the quantum jump 
operators cg^> , acting on pairs of neighboring lattice sites, 
are given by 

c w = (e-^^ + e-^^)^^-^^) 

= (b\ + e-**"' b\, ) (b £ - e**"' b v ) . (4) 

with = cj)£t — (j)£. We specialize to the case that the 
phase difference along the A = 1 . . . d primitive directions 
eA of the lattice is a constant = <p£+e x ~ fyt- The 
quantum jump operators are composed of an annihila- 
tion part, which destroys a particular superposition of 
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FIG. 1: The main panel shows the q = entry of the single- 
particle correlation matrix in momentum space as the inter- 
action strength is varied, obtained from the solution of the 
exact EOMs Q for N = 2, L = 12, J = 0.0, and = 0, 
at the time TfK — 20.0. The inset shows the entire diagonal 
of the correlation matrix for U/k = 0.0 (shaded bars) and 
U/k = 20.0 (empty bars). 



FIG. 2: The main panel shows the entropy in the steady 
state of the system, from the solution of the exact EOMs |l]) 
for the same parameters as Fig. ^ The inset shows the same 
data as the main panel (crosses) logarithmically rescaled with 
respect to the asymptotic value Soo — 0.72 and compared to 
the power-law Soo — S oc U~ a for a < 1.0 (solid line). 



the bosonic wavefunction, and a creation part, which re- 
cycles the bosons after interacting with the bath. While 
the choice of the phase in the annihilation part is crucial 
to the asymptotic behavior of the system (see below), the 
phase in the creation part is a matter of convention and 
could also be omitted. One possible implementation of 
these jump operators in optical superlattices is discussed 
in Ref. [15] . The key point of such choice for the jump op- 
erators is that their nullspace consists of a unique element 
\D) for a fixed number of particles N. It follows that \D) 
is the unique dark state for the whole Liouville operator 
and that the system is attracted to this state indepen- 
dently of the initial conditions, if no competition arises 
from the unitary part of the evolution. In Refs. [15] [23] 
it was shown that the state 

|BEC q ) = (TV!)- 1 ^ ^V iq£ &jj = (N\)~^% N \v^c)(5) 

[a Bose-Einstein condensate (BEC) with momentum q = 
— (j)\e\ in lattice units] is the unique dark state of the 
Liouvillian The uniqueness can be easily under- 

stood in momentum space, where the annihilation part 
of car reads ^a{1 — exp [—i(q\ + (t>\)]}bq. Since the cre- 
ation part of Cij* does not have any eigenvalue on the 
Hilbert space of N — 1 particles, &w features a unique 
dissipative zero mode at the given momentum. The ex- 
istence of a dark state depends here on the jump oper- 
ators being not hermitian. It is important to remark, 
anyway, that the Liouville operator conserves the to- 
tal particle number N = as the jump operators 
just redistribute a particle's superposition and thus ful- 
fill exp (+i(f)N)c££> exp (—icf)N) = ca> . However, the non- 
hermit ian nature of the jump operators entails that the 
dissipative process represented by the Liouvillian ([3| does 



not respect the principle of detailed balance, because the 
probability that the system is projected out of the dark 
state is zero. The dissipative dynamics produced by the 
engineered Liouvillian is necessarily very different from 
a thermalization process and hence the phase portrait 
of the system cannot be understood in terms of a modi- 
fied or perturbed equilibrium picture, but requires a gen- 
uinely new approach. 

An alternative physical picture of the dynamics in the 
presence of the dark state is the following. On two 
neighboring sites along the A direction, the component 
of the wave function that is not in the nullspace of 
bi — e L€ ^ x biJ rGx is mapped onto the whole Hilbert space, in- 
cluding the nullspace, while the nullspace is mapped onto 
itself. Asymptotically, the weight of the wave function 
outside the nullspace vanishes and the phase of the wave 
function between the neighboring sites is locked to 
Since the phase locking takes place on each pair of sites, 
the many-body state is asymptotically projected into the 
BEC state, i.e. a state with long-range phase coherence. 
The dissipative evolution purifies the initial state of the 
system and, in this sense, pumping into the BEC bears 
strong analogies to laser cooling of single atoms [24 . The 
analogy to optical pumping is quite tight, because the 
implementation of the dissipative process involves a co- 
herent (laser) driving element, such that the system is 
indeed driven-dissipative. 



B. Competition between unitary and dissipative 
dynamics 

We start out by considering the effect of strong in- 
teractions on the steady state of the system. In the 
equilibrium phase diagram of the Bose-Einstein Hamil- 
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tonian Q, the quantum phase transition between the 
Mott insulator and the supernuid is understood in terms 
of the competition between the kinetic term proportional 
to J, that tends to delocalize the bosons to decrease the 
energy, and the interaction term proportional to U, that 
is responsible for an increase of the energy when parti- 
cles tunnel between sites and repel each other. To un- 
derstand the interplay of the energy scale in our present 
setup, a novel point of view is required, because the time- 
evolution of the system does not minimize the energy. 
The absence of a variational principle to characterize the 
steady state of the system is indeed a noteworthy dif- 
ference with respect to the equilibrium theory. Still, we 
identify the existence of a competition between the en- 
ergy scales of the interaction U and the dissipation n. To 
this end, we proceed by adopting a rotating frame that 
eliminates the interaction from the Hamiltonian. The lo- 
cal operator that defines the new frame on each site is 
Vi = ex.-p[iUhi(h£ — l)t]. The annihilation operator in 
the new frame reads 

VbtV -1 = exp [— iUht]bi = exp [inUt] \n)i(n\bi . (6) 

n 

The effect of the interaction is to rotate the phase of each 
Fock state differently, destroying the phase coherence of 
the order parameter. Hence we argue that increasing 
the interaction leads to a substantial depletion of the 
condensate and eventually to the disappearance of the 
long-range order in the system. 

Before developing a suitable mean field approximation 
that allows us to tackle the question of the behavior 
of the system in the thermodynamic limit, we provide 
signatures of the competition between interaction and 
dissipation from exact numerical evaluation in a small 
system. For this purpose, we integrate numerically the 
equations of motion (EOMs) ([!]), in the Hilbert space 
for N = 2 particles in L = 12 lattice sites with periodic 
boundary conditions, using a fixed- stepsize fourth-order 
Runge-Kutta method. We use the final density matrix 
p(Tf) to compute the single-particle correlation matrix 
(bqbq') and the entropy S = — Tr[plnp], where q,q' are 
momenta in the discrete Brillouin zone of the finite lat- 
tice. Fig. [I] shows that, as the interaction strength is 
increased, the occupation of the zero-momentum state 
decreases monotonically. This is a signature, in the fi- 
nite system, of the depletion of the condensate fraction 
that we expect in the thermodynamic limit. We also 
see that the distribution of the momenta is localized in 
the absence of interactions and broadens substantially for 
U > ft, analogously to the broadening of the condensate 
peak produced by a finite temperature. Such broaden- 
ing is not due to a rearrangement of the single-particle 
states, e.g. the momentum states being not the appro- 
priate eigenbasis for the system, because, as we show in 
Fig. [2j the entropy steadily increases with the interac- 
tion strength. This means that, in the presence of in- 
teractions, the Liouvillian is no longer able to drive the 
system towards the pure dark state. It is interesting to 



note that the entropy converges towards a definite con- 
stant value as the interaction strength grows, as is clear 
from the analysis in the inset of Fig. [2] This fact suggests 
that the crossover exhibited by the finite system may in 
fact turn into a transition from the condensed state to a 
well-defined new phase. In the following section we val- 
idate this picture and provide a detailed analysis of the 
phase diagram. 

III. NONEQUILIBRIUM PHASE DIAGRAM IN 
THE MEAN FIELD APPROXIMATION 

In order to capture the physics of the above depicted 
nonequilibrium phase transition on a semi-quantitative 
level, we develop a formalism which allows us to describe 
the physics of both the well-controlled limits of weak 
and strong interactions. In the weakly interacting limit 
U/k <C 1 we must recover the condensate phase with 
long-range order found in Ref. [15], while in the strong 
coupling limit U / n ^> 1 we expect a disordered phase 
with diagonal density matrix. To this end, we tackle the 
solution of Eq. ([!]) within a mean field approximation, 
in the form of a Gutzwiller-like ansatz where the system 
density matrix is factorized in position space 

P = <g)Pi- (7) 

t 

The reduced density matrix pi = Tr^ep on the site £ 
is obtained by tracing out all the the other degrees of 
freedom but those of the £th site. The equation of motion 
for the reduced density matrix then reads 

d t pi(t) = -iTr#[H,p(t)]+Tr#C\p(t)] . (8) 

The computation of the traces is reported in App.|A] The 
Hamiltonian part results in a local commutator 

- iTr#[H, p(t)] = -i[h £ (t), Pi (t)] , (9) 

with the effective reduced local Hamiltonian given by 
he = -JYs(ineX$i')b\ + fte'fii) + \Ufn(fn - 1) - phi 
[cf. also Eq. ( |A6[ )]. This is precisely the form obtained in 
the ordinary Gutzwiller approximation for a mean field 
ground state calculation on the Bose-Hubbard model. 
Note, however, that due to the restriction of the ground 
state, the state must be pure and thus the onsite density 
matrices in the latter case are of the form pi = 
In our case we do not have an argument in favor of the 
purity of the steady state for strong interactions so we 
allow for mixed steady states as well. 

The trace of the Liouvillian part, assuming an axial 
geometry in the system (see App.[A|, produces the struc- 
ture 

Tr^£[p(i)] = K J2[^ Pi Ap -ApA r ePi - Pi ApA r e } 

r,s 

x [r^u+i + r& = -i + (z- 2)r& =0 ]U0) 
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where z is the coordination number of the lattice, the the number operator. We reiterate that the treatment of 



vectors of operators A are given in E q. (All), and the 
matrix of coefficients Tg j(J in Eq. (A13). 



At this point we discuss the validity of the factoriza- 
tion approximation ^ in real space. First we note that 
the dimensionality of the system enters the EOMs only 
through the coordination number z of the lattice. In gen- 
eral, one expects higher reliability of mean field approxi- 
mations in higher spatial dimensions. For the (fermionic) 
Hubbard model, arguments based on the central limit 
theorem for spatial dimension d — >• oo lead to the ob- 
servation that mean field theory becomes exact in this 
limit [25 . The statistical nature of the argument sug- 
gest that a similar reasoning could be applied also here. 
The final results of the analysis do not depend qualita- 
tively on the precise value of z hence, when the choice 
of a definite dimensionality is necessary, we consider for 
simplicity a one-dimensional (d = 1) system. Moreover 
we observe that the system is represented only locally as 
a quantum state, i.e. a vector or a density matrix in the 
local bosonic Fock state of each site £. Any correlation 
function (OgOg>) that involves the evaluation of opera- 
tors on different sites simplifies to the product (Og)(Ogf). 
The connected part of the correlation function between 
different sites is entirely neglected. However, the onsite 
correlations are evaluated exactly and thus this method 
is quite suitable to extract informations on a wide range 
of variation of U / k, from the weak to the strong coupling 
limit. In the limit of weak coupling, additional arguments 
in favor of our approximation scheme are available (see 
Sec.lVA). 



Technically, the Gutzwiller approximation transforms 
the full many-body system into a collection of coupled 
single-site systems, that are amenable of a much simpler 
analytical treatment. Moreover, in a numerical solution, 
the total dimension of the state vector of the system de- 
creases from the full size of the Hilbert space to L times 
the size of the local Hilbert space. This makes much 
larger arrays numerically tractable. On the technical 
side, the price to pay for the decrease in size of the Hilbert 
space is that the equation of motion becomes nonlinear 
in the state vector, contrary to the Schrodinger and von 
Neumann equations. More precisely, the coefficients that 
enter the effective local Hamiltonian hg> and the matrix 
r^/ )Cr depend on a set of local correlation functions 



^g = (bg), (if), (hgbg) 



(11) 



evaluated with the density matrix of the sites t neigh- 
bors of i f . The appearance of a nonlinear equation from 
an originally linear equation for the system density op- 
erator in the framework of a mean field approximation 
is actually familiar from e.g. the Gross-Pitaevski mean 
field theory for weakly interacting Bose gases. There, a 
nonlinear classical field equation for the condensate ex- 
pectation value results from a mean field approximation 
on the original linear TV-body Schrodinger equation. Fi- 
nally, the mean field equations conserve the average value 
n of the total particle density but do not commute with 



the correlation functions (11) on each site is fully quan- 
tum mechanical, i.e. there is no approximation on the 
value of the average of any local operator. In particu- 
lar, this procedure yields a theory that contains much 
more information than a Gross-Pitaevski equation [see 
Eq. (13) and the following discussion], where also the 
product of local operators ( pTj ) is factorized, so that all 
available information is encoded in x/j£. 

The effect of the phase <j) in the Liouvillian is easily 
understood for a homogeneous state by considering the 
change of frame of reference bg = e+^bg. In the new 
frame, the Liouvillian maintains the same form, with 
bg H> bg and (f) 0. In the effective Hamiltonian, instead, 
it is necessary to apply bg \-> bg but also J \-> Jcoscj). 
Hence, we can reduce to considering the homogeneous 
states of the system with <\> = only, as the general 
case amounts to a mere reduction of J if <p < tt/2. For 
(j) > 7r/2, however, the dark state of the Liouvillian cor- 
responds to a boosted condensate that is known to be 
dynamically unstable [26] and consequently does not fea- 
ture a stable homogeneous steady state. 



A. Phases of the homogeneous steady state. 

In this section we study the steady states of Eq. (|8| 
and validate the picture, suggested in Sec. |IIB[ of an 
interaction-driven nonequilibrium phase transition that 
destroys the condensation in the system. We first iden- 
tify analytically two exact steady states, for vanishing 
and large interactions, and then provide numerical re- 
sults that interpolate between these two cases. 

Let us consider first the pure coherent state [27] 



oo / v 

v=0 



(12) 



(where v is the index in the local Fock space) and the 
corresponding density matrix = \ip)g{ip\. The Liou- 
villian ( [To] ) vanishes exactly when applied to the homoge- 
neous coherent density matrix, that provides hence a con- 
sistent dark state solution for the dissipative dynamics 
and pinpoints an exact steady state solution of the model 
at J = U = 0. Moreover, for this density matrix, all on- 
site correlation functions factorize in Fock space. As a 
consequence, all the information of the theory is encoded 
in the simplest one, namely the one-point function (bg), 
which is precisely the order parameter ( [TT] ) of the con- 
densate. The equation of motion d t (bg) = Ti[bgd t pg(t)] 
for the order parameter reads explicitly 

- id t ifjg(t) = J ^2 ipg' + pipg - U\^g\ 2 ^g 

+2ik Ml - ^ + - 1^' | 2 V#3) 

(£>\£) 



FIG. 3: Time evolution of the condensate fraction |^| 2 /n on 
the sites 1 < £ < L — 22 of the lattice, according to the mean 
field EOM ([8}, for n = 1.0, J = 0.0, U = 3.0, and = 0. The 
color of the surface represents the phase difference Si between 
neighboring sites (white for 8g = 0, black for 8g — ±7r). 



This equation constitutes a dissipative Gross-Pitaevski 
equation (GPE) and, indeed, the unitary part is identi- 
cal to the standard lattice GPE. If we choose a homo- 
geneous state tp£ = i\) [analogous to the BEC state ([5|], 
the dissipative terms vanish and the equation reduces 
to —id t ^ = (zJ — Uip*i/j + /i)^, i.e. the system under- 
goes a global phase rotation that can be removed with a 
precise choice of the chemical potential \i — —zJ + Un. 
More generally, we will constantly take advantage of the 
possibility to fix the chemical potential and remove resid- 
ual phase rotations, ensuring a time-independent steady 
state as illustrated in this simple case. The role of the 
chemical potential in our nonequilibrium analysis is fur- 
ther discussed below in Sec. IIIIBI The coherent state is 
not a general solution of Eq. (|8| precisely because a differ- 
ent phase frequency appears for each different correlation 
functions but only one can be removed with the chem- 
ical potential. Anyway, the similarity of the Eq. (13) 
to the GPE suggests that small values of the interac- 
tion strength U do not change qualitatively the dynam- 
ics of the system, which is the result rigorously derived 
in Ref. [15] within the Bogoliubov approximation. 

Let us consider now the opposite limit of a diagonal 
density matrix state, in which the coherences are totally 
absent. According to the heuristic argument presented 
in Sec. II B this situation is expected for large values of 
the interaction strength. In this case, the kinetic term of 
the Hamiltonian drops out, since no off-diagonal order is 
present and thus (bi) = 0. The interaction term drops 
out as well, because the operator fig commutes with any 
state that is diagonal in the Fock basis. The only remain- 
ing energy scale in the system is the dissipation strength 
k. Finally, the matrix of coefficients in Eq. (10) reads 
r^ )Cr = diag(0, (n^), (fit) + 1, 1) and the EOM for a homo- 
geneous system reads 



dtpi(i) = 2ft(n + l)(2bepeb\ - b\bepe - pib\bg) 

+2K,n(2b\peb£ - bS\pi ~ pihb\) . (14) 
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FIG. 4: Characterization of the homogeneous steady state 
for J = 1.5ft and = 0. (a) Stroboscopic plot of the time 
evolution of the condensate fraction for n — 1.0, starting from 
an initially fully condensed state, (b) Comparison of the full 
numerical solution of the homogeneous system for n — 0.1 to 
the analytical low density result, (c) Critical point in the ho- 
mogeneous phase diagram as a function of density for different 
values of J = 0.0, 0.5, 1.0, 1.5. 



The latter equation of motion describes precisely a 
bosonic mode bg coupled to a bath with thermal occu- 
pation n = (hi) [27 . It is intriguing that in this case 
n is a property of the system itself and hence we may 
say that the system acts as its own bath. In particular, 
this effective thermal bath, which drives the dynamics, 
is not related to the physical bath that is necessary to 
implement the jump operators Q. In the steady state 
d t p i — the remaining energy (rate) scale n drops out, 
and the only physical scale is the average particle num- 
ber n. The steady state solution is readily seen to be a 
thermal-like state 



4 th) = Ei 

u=0 



(n +l)("+i) 



(15) 



We emphasize that this state emerges despite the fact 
that the system is driven out of thermodynamic equilib- 
rium, not as a consequence of the latter. Finally, we note 
that the nonlinearity of the physically motivated mean 
field approximation allows us to describe two distinct 
dynamical fixed points, namely one ordered state with 
finite off-diagonal expectation values such as (hi), and 
one disorder phase with vanishing coherences, in some 
analogy to 4 -theory for weakly interacting boson sys- 
tems. We also emphasize that the breaking of the U(l) 
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phase symmetry (associated to the conservation of par- 
ticle number) takes place spontaneously in the ordered 
phase, since both Hamiltonian and Liouville operator in 
the original equations ([l}j3| conserve the number of par- 
ticles. 

We turn now to the numerical integration of the 
EOM (|8|, specified to a lattice with L sites and peri- 
odic boundary conditions. The integration is performed 
with a fourth-order Runge-Kutta algorithm, using a fixed 
stepsize typically taken St = 10 -3 ^ -1 [28]. The local 
Hilbert space is truncated to a maximum particle num- 
ber ^ max i$ 20, substantially larger than the typical par- 
ticle densities n < 1.0 that we consider here. Due to the 
periodicity of the lattice, the momenta are confined to a 
discretized Brillouin zone such that the phase <\> assumes 
the values 27rm/L, with m G N. For each site £ in the 
lattice, the execution of the integration returns the local 
correlation functions. To measure the degree of coher- 
ence along the array, we compute the phase differences 
Si between the first-order correlation functions on neigh- 
boring sites 



2irSi/L — Imln 



KWm)l 2 



(16) 



A state with definite momentum is such that Sg is integer 
and constant along the array. A typical output of the 
numerical integration is displayed in Fig. [3J where we 
show the condensate fraction \tpg\ 2 /n as a function of 
time on a lattice with L = 22. The array is initialized as a 
collection of thermal states (15), slightly perturbed with 
tiny random off-diagonal fluctuations. The figure shows 
that the system is driven out of its initial thermal state, 
develops a uniform macroscopic condensate fraction < 
1.0 and, at the same time, the fluctuations of the local 
correlations vanish and a definite phase difference Sg = 
is established throughout the array. The remarkable 
property of the Liouvillian (|3|, namely that local phase 
synchronization results in the establishing of long-range 
order [15 , is visually clear in this picture. 

We use now the numerical integration of the EOM ([sj) 
to compute the asymptotic value of the order parameter 
as the interaction strength is increased. Since we inves- 
tigate the properties of the homogeneous steady state, it 
is sufficient to consider a single site (L = 1). We ver- 
ify that at the end of the time evolution the infinity- 
norm max{|<9tP£ ; ^' | 1 < v, v' < ^ max } converges towards 
zero. In panel (a) of Fig. [4] we plot the condensate ex- 
pectation value as a function of increasing interaction 
strength for different, equally spaced times, starting the 
evolution from an initial fully condensed state (U = 0). 
We show that the asymptotic values of the order param- 
eter interpolate continuously between the expected val- 
ues \ip\ 2 /n = 1 for U = (coherent state) and ip = 
(thermal state). Similarly, the asymptotic values of the 
entropy Si = — Tr[p^lnp^] shown in Fig. [5] interpolate 
between S = for the pure state at U = and the value 
S = 2 In 2 that corresponds to a thermal state with uni- 
tary average density. Between the two extreme values, 



T 



bjO 
O 




-0.5 



0.5 

log 10 U/k 



FIG. 5: Asymptotic value of the entropy S as the interaction 
strength U is increased, for J = 0.0, 4> = 0. The dashed line in 
the main panel corresponds to S ~ (U/k) 1 ' 7 ; the dot-dashed 
line corresponds to S = 2 In 2 ~ 1.38. The final time of the 
time-evolution is Tf k ~ 40.0 in the main panel and Tf k ~ 2.1, 
2.9, 6.4, 15.0 in the inset, from the bottom to the top curve. 
The arrow marks the critical point U c . 



both quantities present a non-analyticity that develops 
in time, as the system converges better and better to the 
steady state. This effect is clearly a qualitative modifica- 
tion over the smooth crossover that appears in the finite 
systems of Fig. [l]and|2j We may conclude that, in the 
thermodynamic limit described by the Gutzwiller ansatz, 
the competition of phase-enhancing dissipation and de- 
phasing due to the interaction leads to non-analytic be- 
havior, i.e. a nonequilibrium phase transition. 

To reinforce the previous considerations, we analyze 
the convergence of the system to the steady state in the 
proximity of the non-analyticity. We verify in Fig.[6]that, 
in this case, the convergence to zero takes place with a 
power law \ip(t)\ ~ £ -1 / 2 , in contrast to the exponen- 
tial convergence generically granted by the dissipative 
dynamics. To be more precise, in the Gutzwiller ap- 
proach spatial correlations are neglected by construction 
and we correspondingly observe exponentially decaying 
correlations, except for the critical point [29J. To in- 
troduce a relation between the time-evolution and the 
criticality of the system, we note that time sets an en- 
ergy scale 1/t for the system, which may be viewed as an 
irrelevant coupling, corresponding to an attractive direc- 
tion of the fixed point of the system tuned to criticality 
by choice of the relevant coupling U/k. Indeed, we ob- 
serve that the non-analyticity corresponding to the crit- 
ical point develops only in the limit where this scale is 
removed. In the vicinity of the critical point the conden- 
sate amplitude scales as ~ exp (— A max t)/£ a , where 
Amax(^//^, J/k) is the real part of the largest eigenvalue 
of the system, vanishing at criticality. There, the system 
shows polynomial behavior with universal critical expo- 
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FIG. 6: Time-evolution of the absolute value of the order 
parameter in a homogeneous system (L = 1) for J = 0.0, 
U/k ~ 7.2, n — 1.0 in linear (main panel) and logarith- 
mic (inset) scale. The dashed line shows a power-law decay 

\*P(t)\ OC 1/y/t. 



FIG. 7: Nonequilibrium phase diagram, featuring thermal- 
like and condensed homogeneous steady states and a domain 
where the condensed state is unstable. The black (blue) line 
is the results of the numerical techniques explained in App.[B] 
and App. [C] for n = 1.0 (n = 0.1). The red line is the ana- 
lytical result in the limit n C 1 and is repeated in the inset 
with a different choice of the coordinates for the axes. 



nent, defined as 



lim — 

t^oo alogl/t 



t^oo OLOgt 



(17) 



which evaluates to a = 1/2 as deduced from Fig. [6j in- 
dependently of J/k or the mean filling n. This is the 
behavior expected within mean field theory and it does 
not depend on the dimensionality of the system. While 
the mean field critical behavior is another clear signa- 
ture of a second order phase transition taking place in 
our nonequilibrium system, it is not sufficient to charac- 
terize the true universality class of the system. For this 
purpose, the long wavelength spatial fluctuations would 
have to be included in an appropriate way, which is not 
possible within a site decoupling Gutzwiller approach. 

Finally, we provide in Fig. [7] the position of the criti- 
cal point as the hopping amplitude J and the interaction 
strength U are varied, for some values of the density n. 
To draw the phase border with good accuracy, it is rather 
awkward to integrate Eq. (|8| repeatedly for long times. 
We resort instead to the analysis of the EOMs linearized 
around the homogeneous thermal state. The key obser- 
vations here is that the thermal-like state (15), contrary 



to the coherent state (12), is always an exact steady state 
of Eq. Its stability property, of course, are different 
in the two phases and this must appear in the spectrum of 
the linearized equations. In the thermal phase, the ther- 
mal state is an attractive fixed point, while it is at least 
a saddle point in the condensed state. Around the stable 
fixed point our approach analyzes the linear response of 
the system to small perturbations, while around a generic 
fixed point it assesses the existence of unstable directions 
in the phase space. In generic nonlinear equations, two or 
more (meta-)stable steady states could occur. This does 
not seem to be the case in the present model based on our 
numerical findings, and is also reasonable given the fact 



that the evolution of the density matrix elements dtpi-vv' 
is generated by a quadratic form. (One could also con- 
ceive a steady state that is unstable, but still is stable at 
any order in the Taylor series: this again does not happen 
here.) The derivation of the linearized EOMs is deferred 
to App. [Bj but we mention here that this method allows 
us to test the stability of several hundred pixels in the 
phase diagram, in a short computational time, truncat- 
ing the Hilbert space to a maximum number of particles 
^max larger than 10 2 . The latter condition is particularly 
important because the diagonal elements of the thermal 
state decrease only as a power law with the Fock index 
v. 

Quite nicely, this approach allows us also to obtain 
the analytical form of the phase boundary in the limit 
n <C 1 of small particle fillin g. T o this end, we con- 



sider the linearized equations (B9) around the thermal 



state with onsite Fock spaces truncated at occupation 
number ^ max = 2, appropriate for sufficiently low filling 
but still capable of describing the effect of interactions. 
Then we have a 3 x 3 linear system where one eigenvalue is 
iU '/2k and corresponds to the v = state. The other two 
eigenvalues A_, A + can be computed straightforwardly. 
Keeping only the first order in the average density, one 
eigenvalue A_ has negative real part while for the other 
Re[A+] = 4n^(4J 2 -/7 2 +4J/7+32^ 2 )/[(2J+/7) 2 +64^ 2 ]. 
The curve in the (J, U) plane where A + = defines the 
phase border, i.e. where the thermal state changes from 
stable to unstable. In the current approximation it reads 



U = 2J+yfKP 



(18) 



with U c = V32ft the critical point at J = (z = 1). This 
expression is shown as a red solid line in the main panel 
and in the inset of Fig. [7[ Below in Sec. IV we point out 
an alternative way of deriving this result. By reducing 
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the analysis of the equations of motion to the first three 
Fock states it is also possible to extrapolate the position 
of the critical point for variable density, as shown in the 
panel (c) of Fig. [4] 

Two features of the phase diagram of Fig. [7] are dra- 
matic manifestations of the system being out of equilib- 
rium. First, even in the limit of vanishing dissipation the 
system presents a transition from a thermal-like to a con- 
densed steady state at J/U = (y/2 - l)/2 ~ 0.21. This 
picture does not connect to the phase transition between 
the Mott state and the superfluid state that character- 
izes the phase diagram of the equilibrium Bose-Hubbard 
system. Formally, such difference is possible because the 
limits k — » and t — )• oo do not commute and defeats 
the naive expectation that a small dissipation applied to 
a system does not wash out entirely its zero-temperature 
phase diagram. A more physical explanation is the fact 
that the engineered dissipation in the Liouvillian Q does 
not respect the principle of detailed balance, as discussed 
above, and hence its action is radically different from a 
bath that induces thermalization. A notable example of 
the latter case is a Josephson junction array coupled to 
a dissipative bath [21 , which suitably chosen can stabi- 
lize the superconducting ordered phase. The profound 
difference with respect to our setup stems from the fact 
that there the system plus bath is in global thermody- 
namic equilibrium and the associated phase transitions 
(thermal or quantum) are equilibrium ones. 

Second, we note the absence of commensurability ef- 
fects. At the mean field level and at zero temperature, 
one expects that any mechanism suppressing superflu- 
idity leads to a Mott state for commensurate filling. In 
fact, given the density matrix p£ = ^ v \ v)t(v\ p v without 
off-diagonal order (which represents long-range order in 
the mean field approximation), the only solution for the 
steady state is a Mott state p v = S u ^ no with quantized 
particle number no, because the conditions ^ v p v = 1 
(normalization) and ^2 u p1 = 1 (purity at zero tempera- 
ture) must hold simultaneously on each site of the lattice. 
In contrast, in our case the dephasing of correlations in- 
duced by the interaction (see discussion in Sec. II B ) leads 
to a diagonal steady state that does not obey an addi- 
tional purity constraint. The dephasing of the coherent 
initial state makes the density matrix more and more di- 
agonal, but does not lead to an additional localization 
in Fock space [30]. We also observe directly that the 
phase diagram in Fig. [7] does not change qualitatively 
from n < 1.0 to n = 1.0. In conclusion, the nonequi- 
librium quantum phase transition presented here shares 
features of equilibrium quantum phase transitions in that 
it is interaction driven, and of classical phase transitions 
in that the ordered phase terminates in a thermal state. 



B. The role of the chemical potential 

In our formulation of the many-body nonequilibrium 
problem, we make use of a chemical potential, which at 



first sight is a concept tightly bound to thermodynamic 
equilibrium. In order to better understand its meaning 
and practical use in the present setting, we would like to 
compare our situation to the thermodynamic equilibrium 
of an interacting Bose gas, say in the continuum, in which 
the temperature is low enough that spontaneous symme- 
try breaking takes place. Here, the role of the chemical 
potential is twofold: first, it conceptually fixes the (aver- 
age) particle number; second, its actual value is chosen 
to fulfill an equilibrium condition, namely to remove ho- 
mogeneous terms which are linear in the fluctuation op- 
erators which occur upon expansion about a condensed 
ground state, a condition which in turn is equivalent to 
the existence of a gapless Goldstone mode. 

In our nonequilibrium case, there is no coincidence of 
these two roles of the chemical potential; its choice is 
only connected to the second one. Indeed, the average 
particle number is an exactly conserved quantity of our 
mean field master equation. Hence, the average particle 
number is fixed via the initial state. However, the intro- 
duction of a chemical potential term is needed to ensure 
the existence of a steady state, in which all time- local cor- 
relation functions become time independent. Mathemat- 
ically, this choice can be made by requiring that the equa- 
tion of motion for the spatially homogeneous fluctuation 
ipi(t) — -0(0) experiences no "driving" for all £, i.e. there 
is no constant term generating its evolution. This con- 
dition on the equation of motion is precisely analogous 
to the condition of the vanishing of the coefficient of the 
term linear in the fluctuation in a Hamiltonian. 



C. Dynamical instability in the nonequilibrium 
phase diagram 

Having classified the homogeneous steady states of the 
system and presented the phase border in Fig. [7| we con- 
sider now the stability of the condensed phase with re- 
spect to inhomogeneous perturbations. Although in gen- 
eral we expect that an initial inhomogeneous system is 
driven towards the homogeneous dark state defined by 
the Liouvillian (as discussed in Sec. [H] and shown for a 
paradigmatic case in Fig. [3| we cannot discard the pos- 
sibility that the competition between unitary and dissi- 
pative dynamics destabilizes the convergence process in 
some range of parameters. The stability of the homoge- 
neous steady state can be tested numerically by subject- 
ing the steady state to a small perturbation and moni- 
toring the dynamics following the perturbation. Due to 
the nonlinearity of the dynamics described by Eq. ([8]), 
the results of this analysis depend on the kind of pertur- 
bation that is applied to the system. From the theory of 
the linear response of many-body systems, we know that 
the spectrum of the collective excitations is obtained by 
studying the density-density response function, hence we 
perturb the density matrix of the system with the addi- 
tion of a coherent component dpi that oscillates periodi- 
cally in space ex cos (2tt£/ L) with the maximum wave- 
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FIG. 8: Time-evolution of the density wave generated by the 
dynamical instability, according to the mean field EOM (g}, 
for L = 22, n = 1.0, J 0.0, and U/k 3.0. The main 
panel shows the fluctuations 5ne of the local density ng with 
respect to the average value (short dash, long dash, solid, from 
the earlier to the later time). The inset shows the absolute 
values of the fluctuations in logarithmic scale (solid) and the 
prediction of the growth rate obtained with the numerical 
linearization of the equations of motion around the steady 
state (dashed). 



length allowed by the size of the lattice. In Fig. 
we show that, for specific choices of the parameters, the 
magnitude of the perturbation 5ri£ is not damped, but in- 
creases exponentially in time as e 7 ^, with 7 > 0, preserv- 
ing its shape. The robustness of the spatial profile Sn(t) 
and the fact that it increases in time with a well-defined 
instability exponent 7 classifies this phenomenon as an 
unstable mode. The system then manifests a dynamical 
instability, characterized by the wavelength of the charge 
density wave show in Fig. [8] A more thorough analytical 
understanding of the origin of this instability is deferred 
to Sec. [V| after a suitable approximation scheme is de- 
veloped. 

To delimitate systematically the range of parameters 
where the dynamical instability takes place, we linearize 
the EOM ([8| around the steady state and study the spec- 
trum of eigenvalues, analogously to the procedure per- 
formed in the previous section to find the border of sta- 
bility for the thermal state. In this case, however, the 
substantial analytical simplifications deduced in the pre- 
vious case cannot be applied because the form of the 
condensed steady state is not known in general and, in 
any case, the EOMs for the entries of the density matrix 
have a much more intricate structure of couplings. Hence 
we resort to the numerical procedure outlined in App.[Cj 
which provides the shape of the unstable domains shown 
in Fig. [7] with dashed lines. We see that an unstable do- 
main exists, that includes the whole J = axis where 
the thermal state is not stable, and is larger as the in- 
teraction strength increases. The phase transition from 
the thermal to the condensed state is then substituted 
by a transition from the thermal to the unstable phase 



FIG. 9: The main panel compares the time-evolution of the 
condensate fraction according to the mean field EOM (|8| and 
the EOM for small density (solid and dashed lines, respec- 
tively) for n = 0.1, 0.03, and 0.01 from bottom to top in a 
homogeneous system (L = 1). The inset shows in a bilog- 
arithmic scale the difference Ac between the results at the 
final time Tf = 20.0 for the three values of the density. 



for small enough values of J. The border between the 
condensed and the unstable regime in the proximity of 
the origin is characterized by a fixed ratio J/U, as can be 
clearly seen in the inset of Fig. [7J To interpret this fact, 
we note that the condensed eigenstate of the Liouv illian 
(10) is also an eigenstate of the kinetic term in Eq. (A6), 



so J is the strength of an unitary term "compatible" with 
the engineered dissipation. On the contrary, interaction 
U and dissipation compete strongly, as discussed in detail 
above. We may argue then that, in general, the steady 
state of the Liouvillian can be made more robust by in- 
creasing the compatible energy scales with respect to the 
competing energy scales in the Hamiltonian. 



IV. ASYMPTOTIC BEHAVIOR AND 
CRITIC ALITY IN THE LIMIT OF SMALL 
DENSITY 

We present now several analytical arguments that 
confirm the nonequilibrium phase diagram presented in 
Fig. [7J We consider the equations of motion for the most 
general normal ordered correlation functions of our sys- 



tem c n 



- /it 



(bgbg). In general, these correlation func- 
tions form an infinite nonlinear coupled hierarchy. How- 
ever, we find that the a priori infinite hierarchy of mean 
field equations for the onsite correlation functions ex- 
hibits a finite subset which decouples from the rest in the 
limit n — >> 0. Importantly, we find that all qualitative ef- 
fects of the nonequilibrium phase diagram are contained 
in this reduced system of equations of motion, making 
this limit an important tool for an analytical understand- 
ing of the physics of our system. In Fig. [9] we show that 
the results obtained with the numerical integration of the 
reduced equations compares quantitatively very well with 
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the full equations for average densities even larger than 
n = 0.1, i.e. a realistic figure for typical fillings in optical 
lattices. 

For the full system of correlation functions, the general 
form for the equation of motion d t c VA ^ in closed form 
is reported in Eq. (Dl). We observe that c PA ^ couples 
to the correlation functions c p+ i 5 i^ and c Pj g+i^, thereby 
not allowing for a truncation to small subset of correla- 
tion functions a priori. However, progress can be made 
by introducing a power counting bg ~ n 1//2 . This gives 
an upper bound on the true powers of the operators, as 
correlation functions in the state without phase symme- 
try breaking with p ^ q vanish exactly. The leading 
power for the time evolution of a correlation function is 
consequently given by c P: q : £ ~ n b+g)/ 2 . Writing down 



the system of equations of motions (11), it is then obvi- 



ous that this system closes if we truncate the equations 
of motion for (6|), (figbi) each of the equations at their 
leading power, thereby physically corresponding to a low 
density expansion. 

We use this result first to obtain analytically the criti- 
cal exponent 7 discussed above. For a homogeneous sys- 
tem with J = the equations of motion read 



d t (b £ ) = mC7(^) + (-iC7 + 4/€)(n^)-4^(M*(6|) 
dt{htbi) — 8nft(6^) + (—iU + inU — 8ft) (n^) 
d t (bj) = (-i/7 + 2m/7-8ft)(^) + 8ft(^) 2 , (19) 

with the chemical potential \i = nU resulting from 
the steady state condition, cf. the discussion following 
Eq. (13) [31 . The exponential convergence of the sys- 
tem to its steady state is described by the linear con- 
tribution, which corresponds to a "mass" or "gap" in 
the jargon of statistical mechanics. We are interested 
in the critical regime of the dynamics, in which the ex- 
ponential convergence is suppressed due to the vanish- 
ing of the gap term and the system is driven to its fixed 



point by the residual nonlinear contributions to the equa- 
tions of motion. We focus on the linear equations for 
(be) and (hebe), with eigenvalues Ai, A2 and we consider 
(b 2 ) as a nonlinear driving term. In the critical regime, 
in which Ai = 0, the secular equation for the eigen- 
value fixes the critical value of the interaction strength 
U 2 — 32ft 2 /(l — n) [that coincides with U c defined after 



Eq. (18) for n <C 1]. The eigenvector corresponding to 
Ai is \JU — 4ft, inU) and hence the EOM for the quantity 
(fi = (iU — 4ft) (b i) + inU(h£bi) does not feature linear 
terms and reads d t ^p\ = — 4ft — 4ft)(^)*(6 2 ). A similar 
quantity if 2 built with the eigenvector corresponding to 
A2, instead, features a regular exponential decay in time. 
The time-evolution of (bg) is given in general by a linear 
combination cpi and y>2 but, for long times, the expo- 
nentially decaying contribution of the latter is negligible 
and hence we can substitute (b \) ~ <fi/(iU — 4ft) into the 
EOM for (fi. At this point we make an adiabatic approx- 
imation by solving the equation dt(b 2 ) = for (b 2 ) and 
substituting the result into the EOM for tpi, that yields 
<9t|^i| 2 = — |(/?i| 4 /(3ft). From the solution \(fi\ 2 = 3ft/t, 
using the expressions above, we finally demonstrate the 
critical behavior \(be)\ = l/(4v^ft). 

We focus now on the asymptotic behavior of the system 
to demonstrate the depletion of the condensate fraction 
shown in Fig. [4] To this end, it is more convenient to 
introduce the fluctuation operators 5b t = bg — ip^ where 
^oo is the (unknown) order parameter in the homoge- 
neous steady state. The EOMs for the "connected" cor- 
relation functions built with the fluctuation operator (see 
App.|D| decouple naturally into linear equations that in- 
volve the connected parts of the mean fields ( pT| , the 
density fluctuation (5b\5be), and as a parameter. In 
the homogeneous steady state, in which (5b i) = (see 
App. [Dl) we have 



d t (5b e ) = -iU^oo(6b 2 e ) - 2(iU - 4ft)^oo(^^) - (iU - An)(5b\5bf) - iU^ + z/ii^ = 

dt(Sb\Sbi) = iU^iSbj) - iU^iSbj)* + (iU - 4ft)^oo(^^ 2 ) - (iU + ^^(J&Jtof )* - lQn^(5b\5k) = 

d t (5bj) = -(iU + 8K + UJ)(5i%)-iUil>Z o =0 

d t (5b\5bj) = -(iU + 8ft + 2iJ)(5b\5bj) - 2(iU + A^^(5b\5ky = . (20) 



The nonlinearity of the theory is now entirely contained 
in the identity n = (5b\5bg)\r l p oo — |^oo| 2 ? whence the value 
of the order parameter is obtained once the other equa- 
tions have been solved. Practically, in the solution we can 
assume that the order parameter is real without loss of 
generality and fix the chemical potential to ensure that 
the fluctuation of the density has vanishing imaginary 



part. The full result for the depletion then reads 



1 



U 2 (U 2 + 4JU + 4 J 2 + 64/t 2 ) 



8[U 2 c 2 + Uci + 8J 4 + 96J 2 k 2 + 128k 4 ] ' 

(21) 

with c 2 = 4k 2 + 3 J 2 and a = 12 J 3 + MJk 2 . For J = 
(z = 1) this reduces to the simple inverted parabolic 
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profile 



J=0 



(22) 



These analytical results are compared to the numerical 
solution of the equations of motion in the panel (b) of 
Fig. |4j demonstrating an excellent agreement between 
the two computations. The quadratic behavior of the 
depletion may be contrasted with the thermal deple- 
tion in a free equilibrium Bose-Einstein condensate in 
three dimensions, where \ip\ 2 /n = 1 — (T/T c ) 3 / 2 with 
T c = 27r/M(n/C(3/2)) 2 / 3 with M the particle's mass and 
£ the Riemann zeta function. The quadratic behavior is 
reminiscent of a harmonically trapped Bose-Einstein in 
two dimensions. It is noteworthy that the phase border 
that one obtains by requiring = in Eq. (21) coin- 
cides exactly with the small-density result in Eq. (18) in 
the linear response strategy of Sec. |III A 




FIG. 10: Real part of the eigenvalues A = — iuo + 7 of the 
linearized EOM for J = 0, n = 0.1, and U = 1.0k. 



V. FLUCTUATIONS-INDUCED DYNAMICAL 
INSTABILITY 

A. Linear response 

We are now in the position of providing a clear- 
cut explanation of the dynamical instability unveiled in 
Sec. |III C] transposing at the level of the correlation func- 
tions the instability analysis around the steady state used 
to delineate the unstable domains in Fig. [7| The start- 
ing point are the EOMs for the connected correlation 
functions outlined in App. [DJ and the knowledge of the 
homogeneous steady state p^ demonstrated in the pre- 
vious section. Then we expand the density matrix of the 
system as pt(t) = p^ + Api(t) and consider the fluctua- 
tions of the connected correlation functions 



(5b\ p 5b\) 



Ti[8b 



£ uu £/00 



A(8b\*5b«) 

]+Tr[Sb\ p Sb q i A Pl {t)] .(23) 



The EOMs are then linearized in the fluctuations 
A(Sbt p Sb ( j) and the values of (5b\ p 5b q ) oc, given by 
Eq. ( |20| ), are such that all the driving (i.e. constant) 
terms vanish upon the correct choice of the chemical po- 
tential outlined above. We remark that the fluctuation 
A(Sbi) does not vanish as the system is not in its steady 
state (see App. |d| and that local fluctuations A(Sb\Sbi) 
of the density modify the flat distribution (hi) = n of 
the steady state. Then we take the Fourier transform 
of the EOMs and rewrite the linear system in terms of 
the connected correlation functions in momentum space 
A(Stf p Sb q ) q oc T,e eiq£ ( s bi Ps i>i)- It is important to no- 
tice that the Fourier transform of the correlation func- 
tions is not simply related to the correlation functions of 
the operators in momentum space, except for the first 
order correlation where it holds A(5b) q = A(b q ) and 
A(5b) q = A(b- q ). (We may denote such correlations 



Aip q and A^_ g , respectively, because the fluctuations of 
the order parameter vanish on the steady state by con- 
struction.) Since the instability shown in Fig. [8] takes 
place at small momenta, in performing the Fourier trans- 
form we focus on the central region of the Brillouin zone 
and substitute the occurrences of the discrete Laplacian 
A^u = U£+i — 2u£ + U£-i with the parabolic dispersion 
-q 2 u q . 

The linear system of EOMs takes the form of a 7 x 7 
matrix (the 3 complex correlation functions, their com- 
plex conjugates, and the real density fluctuation) whose 
eigenvalues A = — iuj+j give the ^-dependent spectrum of 
the system. The eigenvectors of the system correspond to 
modes that evolve as 5pi(t) = (5/^(0)e -zu;£ e +7 *, which are 
stable (unstable) if 7 < (7 > 0). The real part of the 
spectrum for a typical choices of parameters within the 
unstable domain is shown in Fig. [Toj The spectrum fea- 
tures: i) two doubly-degenerate strongly decaying modes 
~ 9.0) that project mainly on the third-order corre- 
lation functions; ii) one decaying mode (7/ft ~ 1.5) that 
projects mainly on the density fluctuation; iii) two low- 
lying modes generated by an admixture of the first-order 
correlation functions. The latter modes are magnified in 



the inset of Fig. 10 The lower mode gives 7 > in a small 
interval around q = 0, hence proving the existence of un- 
stable modes with well-defined momentum. The domain 
where an unstable mode exists in this approximation is 
delimited in Fig. [7] by a red dashed line. 

In general, the decay rate of the modes i) and ii) is 
proportional to O(k) and O(^n), respectively, as it ap- 
pears from inspection of the linearized EOMs. The clear 
separation of the dissipative part of the spectrum into 
groups of modes that have largely different decay rate at 
low momenta suggests that an adiabatic elimination of 
the fastest modes can be performed to bring the 7x7 
linear system in a more compact form. In this way we 
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obtain a renormalized spectrum of the weakly dissipative 
single particle modes, where the instability is encoded. In 
general, the adiabatic elimination in a system 



d t u F = F[u F ,us], d t u s = G[u F ,u s ] 



(24) 



with fast (u F ) and slow (us) modes consists of solv- 
ing F[u F ,us] = for u F and using the result into the 
second equation that becomes dtus = G[F _1 [0, us], us]. 
The application of the procedure introduces new terms 
F _1 [0,^s] that renormalize the equation of the slow 
modes. We apply the procedure once to eliminate the 7 ~ 
O(k) modes and then again to eliminate the 7 ~ O(nn) 
modes. Since the border of the unstable domain extends 
to the origin of the phase diagram in J, U (see Fig. [7]), to 
understand the phenomenon underlying the instability it 
is enough to perform the algebraic manipulations to the 
first order in J and U . We obtain a renormalized 2x2 
linear system for the time derivative dt(A^ q: A^* ) of 
the fluctuations of the order parameter in time, which 
reads 



-i(nU + e q ) 

+inU + s* 



-inU 



+i(nU ■ 



Kq + Tq 



Here, 



(25) 

2nq 2 (2n + 1) is the "bare" quadratic decay 
rate that follows from the analysis of the linear corre- 
lations only, for small interaction and nonzero hopping 
amplitude, and is shown as a black solid line in the 
inset of Fig. [To] e q = Jq 2 is the low momentum ki- 
netic energy. Finally, r q — q 2 (15nJU/ k + 22inU) / 32 and 
s q — —q 2 (nJU/n + 7inU)/16 are the terms that renor- 
malize the slow modes obtained by the adiabatic elimi- 
nation. Without the renormalizing terms the 2x2 sys- 
tem displays the structure of a Bogoliubov equation for 
the condensate modes, with diagonal dissipation n q . We 
point out that a standard quadratic theory can repro- 
duce the Bogoliubov-type EOM but necessarily misses 
the renormalizing terms that are due to third-order lo- 
cal correlations, and thereby the entire physics of the dy- 
namical instability. The latter is thus a clear fluctuation- 
induced beyond- mean- field effect. The eigenvalue of the 
linear system, which approximates the lower mode in the 



inset of Fig. 10, reads explicitly 



-iu q +j q = -iq v / nU(8J - 9nU) /2-n q +lbq 2 nJU '/ '(32/c) , 

(26) 

If J > SnU/9 we can identify c = ^nU(SJ - 9nU)/2 as 
the speed of sound u q = c\q\ of the dissipatively-created 
condensate and we also find a modified decay rate for 
the modes that is quadratic in the momentum. However, 
as J increases, the square root becomes imaginary, the 
contribution of the dispersion to u q vanishes and the de- 
cay rate of the modes is modified by a non-analytic term 
~ |g|, which is positive and dominates over the contribu- 
tion ~ q 2 at low momentum. 

The resulting physical picture for the dynamical in- 
stability in the dissipatively-driven condensate is then 
an interplay of short time quantum and long time clas- 



sical fluctuations. Purely local, i.e. temporal (or quan- 
tum) fluctuations captured in the Gutzwiller approach 
renormalize the complex spectrum at short times. This 
prepares the ground for spatial (or classical) fluctuations 
taking over at later times and responsible for the dynam- 
ical instability. We remark that the Gutzwiller approx- 
imation that we adopt is suited to describe both i) the 
onsite fluctuations that are described without approxi- 
mation in a fully quantum local Hilbert space; ii) long- 
wavelength fluctuations of the coherence that receive the 
main contribution from the disconnected part, which is 
computed solely in terms of the order parameter. 

It is important to realize that the instability arises 
at weak coupling already, where the system is well de- 
scribed by the inhomogeneous Gutzwiller mean field the- 
ory. Without referring to any approximation, it is clear 
that the instability is due to a renormalization of the 
single particle (complex) excitation spectrum, and thus 



encoded in the evolution of (Aip q ,A^* 



The exact 



equation of motion for these quantities is a nonlinear 
equation, with in general nonlocal spatial correlations. 
The Gutzwiller approximation factorizes the correlations 
functions in real space, but treats onsite correlations ex- 
actly. The factorization in real space is justified at weak 
coupling (large condensate) because the dominant scat- 
tering processes are those for opposite momenta off the 
macroscopically occupied condensate. However, the fac- 
torization of the onsite corre lation functions (as in the 
GPE discussed in Sec. Ill A) would not allow to cap- 
ture the physics of the dynamical instability. In sum- 
mary, because the dynamical instability is a weak cou- 
pling phenomenon where our approximations can be rig- 
orously justified, we can rule out the possibility of the 
effect being an artifact of approximations. Of course, the 
above analysis has the status of a linear response theory, 
and thus is not able to shed light the fate of the system 
in the dynamically unstable regime; we will address this 
question in Sec. |VB| 

We note now that the exact Liouville operator ([3| fea- 
tures the dark state also in the case of a single particle on 
the lattice, corresponding to the limit n —> with N = 1 
and L — >> 00. Our approximation yields correctly the ex- 
istence of the dark state at q = 0, since the expression 
for the bare quadratic decay rate reads n q = 2nq 2 in the 
limit n —> 0. On the contrary, the Bogoliubov approxi- 
mation used in Ref. [15] gives K q ~ Annq 2 and hence, in 
the limit n — >• 0, does not yield the correct unique dark 
state. The discrepancy between the two methods lies in 
the fact that, in the Bogoliubov approximation, it holds 
(bS\) — (b\bi) leading to 2n + 1 ~ 2n. The additional 
occurrence of unity that is neglected in the Bogoliubov 
approximation is indeed essential to the existence of the 
dark state and is correctly captured by the Gutzwiller 
approximation. Finally, we remark that, in contrast to 
the dissipative part, the bare theory for the unitary part 
of the time evolution, where two-point correlations are 
not included, coincides with the GPE augmented with 
linear fluctuations. This is rooted in the precise nature 
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of the nonlinearities occurring in the dissipative and uni- 
tary part. 



B. Spontaneous breaking of translational 
symmetry at large times 

Because of the dynamical instability demonstrated 
above, the assumption of a homogeneous steady state is 
not always compatible with the physics of this nonequi- 
librium system. The linear response method applied to 
the onset of the instability suggests that a charge density 
wave (CDW), with well-defined wavelength (see Fig. [8|, 
sets in at later times. To confirm this picture, we proceed 
now to the numerical solution of the nonlinear EOMs in 
the limit of small density. In this case the computational 
complexity of the problem is greatly reduced, as each site 
is described by four numbers only. The EOMs can be in- 
tegrated quickly and reliably on very large lattices with 
L oc O(10 3 ) and for extremely long times tn oc O(10 3 ). 

The spatial profile of the particle density is shown in 



Fig. 11 for a typical time-evolution. We clearly observe 
that a stationary CDW establishes at large times. The 
wavelength Acdw of the periodic modulation is the same 
as the most unstable mode in the lowest-lying branch 
of the dissipative spectrum (the minima in the inset of 
Fig. 10). In other words, the fastest-growing mode in the 



continuum of unstable modes becomes dominant at large 
times [32], its wavelength is transmuted into the wave- 
length of the CDW, and, arguably, the remaining modes 
are eventually suppressed by nonlinear effects that are 
not described by the linear response approach used in 
Sec. |V A| This dynamically stable feature, with a relative 
density modulation (n max — n m i n )/n c± 40% in Fig. 11 
gives an unambiguous, robust experimental signature of 
the fluctuation induced dynamical instability. Typical 
optical lattice setups are smaller than L = 800 sites con- 
sidered here, so it may be difficult in practice to obtain a 
periodic modulation if the most unstable mode has a very 
small momentum. For a lattice with e.g. L ~ 10 2 sites, 
one would find a CDW where only a single wavelength 
fits into the optical lattice, and the measurement of the 
average density on the two halves of the lattice would be 
enough to detect the density fluctuation. The reason to 
focus on larger lattices in our theoretical considerations 
is to exclude that we deal with a finite size effect. Indeed, 
since only a small interval of momenta is dynamically un- 
stable, it may happen for a lattice with 20-100 sites that 
only the lowest momentum q = 2tt/L (in lattice units) 
in the discretized Brillouin zone is unstable (we reiterate 
that the dark state with q = is always stable). So, 
reducing to small lattices only, one cannot rule out that 
Acdw = 1/#cdw °t L, which would lead to a vanishing 
effect in the thermodynamic limit L — >• oo. However, 



Fig. [TT] demonstrates that the spatial scale of the density 
modulation remains finite in the thermodynamic limit 
and coincides with the wavelength of the most unstable 
mode. 
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FIG. 11: Time-evolution of the rescaled density ni(t)/n of 
the system, according to the reduced EOMs for small density, 
for L 800, n 0.1, J 0.0, and U 3.0. The system 
is initially prepared in a fully-condensed state and converges 
to a homogeneous steady state, which becomes unstable at 
k < 2 x 10 3 . A persistent oscillation in the density pro- 
file is visible for later times. The double arrow shows the 
wavelength Acdw — 120 (in lattice units) of the instability as 
computed from the theory in Sec. |V A| 



Intriguingly, the density profile presented in Fig. 11 
suggests that the system in the steady state sponta- 
neously breaks the translational symmetry, which both 
Hamilton and Liouville operator present. This adds to 
the breaking of the global phase symmetry due to the 
presence of the condensate. A state where both symme- 
tries are broken is often termed a "supersolid" although, 
in the equilibrium case, the breaking of the translational 
symmetry is stronger, with a spatial modulation on the 
order of the lattice constant [33] . The spontaneous break- 
ing of the translational symmetry discriminates the dy- 
namical instability considered here from the dynamical 
instability suffered by an ensemble of interacting bosons 
moving with constant momentum through a lattice [26] . 
In the latter case, if the momentum of the atomic en- 
semble is larger than an interaction-dependent thresh- 
old, the atomic current decays and the order in momen- 
tum space is lost. The origin of this instability lies in 
the single-particle dispersion, that becomes negative in 
the outer regions of the Brillouin zone. In our case, in- 
stead, the system occupies a momentum state q = in 
which the single-particle dispersion is positive but the 
long-wavelength fluctuations of a collective (condensate) 
mode are destabilized by local fluctuations. Moreover, 
in the case of Ref. [26] the instability is purely classical 
because it stems from spatial fluctuations. On the con- 
trary, the local fluctuations that play the crucial role in 
our case are genuinely quantum, i.e. temporal fluctua- 
tions. Finally, while the order in momentum space is lost 
in the unstable atomic current, the external drive in our 
case allows the persistence of a (spontaneously) ordered 
pattern, in analogy to synergetic phenomena [34] . 
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FIG. 12: The collapse and revival of the condensate fraction, 
evolved according to the EOM for small density, for L — 32, 
n — 0.01, and J — U — 0.0. The initial state is a collection 
of local coherent states with random phases. The phase is 
changed abruptly from to 3 x 2tv/L at (tft)* = 20.0, after 
the system has reached the steady state. The inset shows the 
probability density V for the time (£k)hw at which the revival 
of the condensate fraction on the first site of the lattice reaches 
0.5 (arrow in the main panel), obtained with 10 3 evolutions 
with random initial phases. 



VI. COLLAPSE AND REVIVAL FOLLOWING A 
PHASE QUENCH 

In the previous sections we have considered the dy- 
namics of the system in the presence of a stationary dis- 
sipative drive. In the spirit of linear response theory, 
this has allowed us to study the dynamics of the sys- 
tem in terms of the instability of a collective mode. In 
this section we focus on the response of the condensate 
to an abrupt change in the conditions of the drive. The 
necessity to adapt to the widely different external con- 
ditions triggers a global response in the condensate and 
we observe intriguing nonlinear dissipative dynamics that 
can be understood in terms of collective variables. Here 
we concentrate on the characteristics of this dynamics in 
the absence of unitary contributions to the time evolu- 
tion, i.e. J = U = 0, with homogeneous density profile. 
We first let the system converge to the steady state cor- 
responding to <j) = 0, i.e. a homogeneous condensate in 
the zero momentum mode, and then we change abruptly 
the phase <fi in the Liouville operator at some finite time. 
In the presence of a finite value of 0, the dark state is a 
condensate with finite lattice momentum —<\> and macro- 
scopic inhomogeneous wavefunction ip£ oc e 1 ^ . We study 
the equilibration dynamics towards the new dark state. 
Our considerations apply to the mesoscopic dynamics [35] 
in periodic lattices, as implemented e.g. in Ref. [36] for 
purely coherent dynamics. The dynamical effects that 
we identify are related to the finite size of these systems. 



The time evolution is reported in Fig. 12 for a typical 
case. The overall dynamics of the system consists of the 
rearrangement of the condensate fraction from the zero 



FIG. 13: The Fourier spectrum \0 q \ 2 of the phases of the 
fluctuations (be) — \ij)i\e ldl , for the same evolution shown in 
Fig. |12| The Fourier spectrum is normalized to its maximum 
at each distinct value of the time. Since the distribution of 
phases 0i is real, a symmetric peak appears in the negative 
part of the q axis and is not shown in the plot. 



momentum of the initial state to the momentum — <j> im- 
printed by the dissipative mechanism after the quench of 
the phase. The rearrangement consists of an initial quick 
collapse of the condensate fraction and a subsequent re- 
vival with different momentum. It is remarkable that the 
revival takes place after a substantial but well-defined 
time interval. Indeed, in the particular case shown in the 



inset of Fig. 12, the fluctuation of the time scale for the 
revival is of the order of 5%. To understand the dynam- 
ical processes that take place during the transient of the 
rearrangement is the goal of the present section. 

Immediately after the quench, the condensate fraction 
decreases exponentially towards negligible values. To de- 
scribe the exponential decay we consider the EOMs (19) 



and we keep only the terms that are linear in the powers 
of bosonic operators, that are dominant in the regime in 
which the correlations become negligible, obtaining 



ft(S*) = -4«[l-cos(0)]&) . 



(27) 



We see that the initial state is depleted with a rate that 
increases as the final phase (j) is quenched further from 
the initial value (j) = 0, where the initial state is stable. 
To describe the revival we consider again the EOMs 



(19) and it is convenient to use the moving frame of ref- 



erence introduced in Sec. |III| The steady state of the 
system is a homogeneous condensate with respect to the 
correlations measured by the transformed operators bg. 
From the analysis of the instability performed in Sec. [V] 
we know that the effect of higher-order correlations is not 
necessarily negligible during the growth of a condensed 
mode. Still, we can neglect the product of correlation 
functions, that are much smaller than the correlations 
themselves before the revival has taken place. Taking 
the Fourier transform of the correlation functions with 
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such prescription gives the coupled linear equations 



1 

2~K 



dt 



(b)g\ _ /2(cos 9 -l) 2cos<A / (b) q 



(nb) q 



4n 



(fib), 



■ (28) 



Computing the eigenvalues of the matrix that defines the 
linear system we see that the growing mode that origi- 
nates the revival can be located only at q 2 < n, that 
in the original frame is a small interval centered around 
— (j). We note that reducing the equations of motion to the 
terms proportional to bg only, as we did to describe the 
collapse of the condensate, is not appropriate in this case 
since the q = — (j) mode would result marginally stable, 
but not unstable. Hence, we identify a dynamical selec- 
tion mechanism [37] that promotes the fluctuations in a 
certain interval of finite width. It would be difficult to 
imagine a scenario in which a single mode is unstable, due 
to the continuity in q of the equations of motion, which 
cannot give discontinuous results at any finite time. The 
modes different from the dark state will be damped sub- 
sequently due to the action of the Liouvillian. 

We turn now to the more interesting question, how 
the mode with lattice momentum —0 is initially popu- 
lated during the time interval in which the condensate is 
depleted. In other words, we are interested in the dynam- 
ics of the condensate fraction ip£ = (6^), which exhibits 
tiny fluctuations around the fully-collapsed value ip£ = 
and which reorganizes at ~ (k)nw into the new reviv- 
ing mode. One naive expectation is that the phases 0£ of 
the coherent fluctuations ipi(t) = \^i{t)\e ie ^ random- 
ize and that the component q with momentum q = — (j) 
increases independently due to the process of dynami- 
cal selection outlined above. On the contrary, it is vis- 



ible from the Fourier spectra shown in Fig. 13 that the 
phases of the fluctuations exhibit a very organized dy- 
namical phenomenon, featuring metastable steady states 
and rapid jumps into new configurations. The Fourier 
spectrum has always a very clear peak that samples all 
the discrete points in the Brillouin zone, from the initial 
q = to the final q = —(j). This contrasts indeed with 
the simpler organization of the phases that takes place 
starting from a random distribution, as can be seen in 
Fig. [3] or for very short times in Fig. 13 



In the following we argue that the momentum drift 
in discrete jumps shown in Fig. [13] can be under- 
stood as a cascade of macroscopic quantum tunneling 
(MQT) events between unstable steady states of the 
time-evolution. Such picture is a qualitative modifica- 
tion of the well-known tunneling from a false vacuum 
into a global ground state that takes place in a system at 
quasi-equilibrium [38]. To justify our picture, we study 
the dynamics of the phases 0£ only, given that the ab- 
solute value of the coherences ip£ is physically negligible. 
The equation of motion for the phases reads, in general, 
id t 0e(t) = (2t/j£)~ 1 d t ^£-(2^)~ 1 d t ^ and we specify this 
expression by keeping only the following particular term 
in the EOM for ip £ 

T^dtMt) = e^+iW - 2Mt) + e-^-i • (29) 



The reason behind this choice is that, to the extent of 
identifying the metastable states of the phase dynamics, 
the linear approximation is enough. Indeed, the linear 
term in the equations of motion is responsible for the 
exponential damping of the variables as we discussed in 



Sec. IV in the context of critical dynamics, or at the 
beginning of the present section where the collapse of 
the condensate fraction is demonstrated. The equation 
of motion for the phases then takes the form 



1 

2k 



d t 0i(t) 



2sin('^?^cos(^ + i 



H-l 



2$) 



(30) 

where is the discrete Laplace operator. The struc- 
ture of this equation is reminiscent of an overdamped 
string displaced by 0£ at the position I. The kinetic term 
oc d 2 is absent from the equation and the first derivative 
represents a friction term. We conjecture that includ- 
ing the two neglected degrees of freedom (i.e. the phases 
of (nibi) and (of)) back into Eq. (29) would reinstate a 
kinetic term into Eq. (30) and couple the phases to an 



effectively fluctuating field, although a more thoroughly 
designed approach (e.g. within the Keldysh action for- 
malism) would be necessary to maintain this assertion. 

Following Ref. [39 , let us parameterize the spatial de- 
pendence of the phases in terms of an "instanton" or 
phase slip configuration, 0t(t) = Q{t)(l — l)/{L — l) + 50£, 
to distinguish the local dynamics 50 £ from the global 
variable Q(t) describing the winding of the phases in 
the circular geometry. It is important to notice that 
the difference 50£+\ — 50 £ between phases on neighbor- 
ing sites is defined in [— 7r/2, tt/2], to minimize the to- 
tal variation of the phase profile. With this convention 
for the difference between phases, the winding number 
reads w = mod[(5,27r]. We first consider the case that 
the global variable is constant and study the dynam- 
ics of the spatial fluctuations, that move according to 
(2K)- 1 d t 50£ ~ cos [(f) + <p] A £ 50 with <p = Q/(L-1). This 
is exactly the equation of an overdamped string, slightly 
displaced from its equilibrium configuration, hence all the 
fluctuations 50 £ converge to zero. We consider now the 
dynamics of the global variable Q(t). If < Q < 7r, due 
to the boundary conditions we have that 0l+i — 0l = 
—Q. The equation of motion for Q(t) reads in this case 
(2f<i)~ 1 dtQ(t) ~ —cos[(j)]Q and hence Q decreases to- 
wards zero (we neglect terms that go to zero for large 
lattices). If 7r < Q < 2tt, the boundary condition 
now is such that 0l+i — 0l = 2tt — Q, that entails 
(2^)- 1 ^[2tt - Q(t)\ ~ -cos[0][2tt - Q(t)} and hence Q 
increases towards 2ir. The same procedure applies to the 
other periodic replicas of the interval [0, 2tt) and shows 
that the global variable Q(t) tends to the closer multiple 
of 2tt. In conclusion, we have shown that the state with Q 
multiple of 2tt (corresponding to a definite momentum in 
the distribution of the phases) is stable both with respect 
to local fluctuation of phases and to a global change of 
momentum. We could hence think of a momentum state 
as a local minimum, whereto the system converges by an 
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FIG. 14: The average phase Oai, for the same evolution shown 
in Fig. |12| At each time, the winding number is given by the 
number of abrupt jumps in the average phase. The locations 
of the phase jumps remain stationary in space, between the 
phase slip events, during which the winding number increases 
by one. 



overdamped motion. 

Having clarified the nature of the metastable configu- 
rations seen in Fig. [l3j we improve our conjecture on the 
effect of the terms that have been neglected in the equa- 
tion of motion for 0£. They necessarily provide a mech- 
anism for the system to tunnel out of such minimum to- 
wards a nearby momentum configuration. The picture of 
a sequence of tunneling events, each with a well-defined 
lifetime, explains the stability of the revival time over 
many configurations (see inset of Fig. 12) and justifies a 
bell-shaped probability distribution for this quantity on 
the basis of the central limit theorem. 

To further clarify the phase dynamics, we follow 
Ref. [39] and compute the average phase difference 6/\£ = 
Arg(6{^ + Af) between sites at distance A£. For 
the initial configuration with zero momentum this quan- 
tity is zero everywhere. We see that the initial phase 
configuration remains stable for a noteworthy time lag 
20.0 < tn < 35.0 following the phase quench at tn ~ 20.0. 
Then a fluctuation in the phases takes place and culmi- 
nates in a phase slip at tn ~ 40.0, that increases the 
winding number of the phase distribution by one. An- 
other phase slip takes place at tn > 40.0 and the final 
one at tn > 120.0, corresponding to the abrupt jumps in 
the peaks of the Fourier spectrum shown in Fig. [13] We 
notice that the dynamics in momentum space takes place 
abruptly while it has a smooth appearance in real space: 
this discrepancy is reminiscent of the instanton descrip- 
tion of the MQT phenomenon, in which the change of 
quantum numbers takes place instantaneously in imagi- 
nary time [39] . 

We point out that the strong suppression of the con- 
densate fraction, which takes place for very small densi- 
ties, is not a necessary condition for the metastable be- 
havior investigated above. In Fig. 15 we report a typical 



FIG. 15: Partial revival of the condensate fraction following a 
phase quench, obtained by integration of the reduced EOMs 
for small density, for L — 32, n — 0.1, and J — U — 0.0. 
The phase cj) is changed abruptly from to 3 x 2tv/L at 
(tK,)* = 20.0. For tK < 50.0 the system stabilizes around the 
momentum 2 x 2tt / L and develops a reduced metastable con- 
densate fraction. Nonetheless, convergence to the steady state 
with momentum 3 x 2tt/L is expected to take place at times 
larger than Tf k = 500.0. The inset shows the condensate frac- 
tion at Tf as the density n is varied. For n < 0.05 the system 
reaches the steady state and the condensate revives fully (see 
also Fig. 12). For 0.05 < n < 0.45 the time-evolution is simi- 
lar to the main panel, and the size of the condensate fraction 
depends strongly on the density. For 0.45 < n the system 
develops a metastable condensate fraction with momentum 
2tt/L. 



time-evolution for larger densities (still well described by 
the reduced EOM for small density, see Fig. |9| in which 
the condensate fraction features a partial revival that 
persists as long as the true steady state is not reached. 
While the explanation of the collapse of the condensate 
still holds, we do not have at present a dynamical model 
for the generation and persistence of a reduced conden- 
sate fraction. As the density in the lattice is increased, 
the momentum of the metastable state changes and, at 
the transition point, a discontinuity in the fraction of 
revived condensate takes place. To reproduce a curve 
like that shown in the inset of Fig. [9] it is only neces- 
sary, in principle, to measure the condensate fraction at 
time tK > 10 2 . Hence, these features provide quite clear 
qualitative observable signatures of the underlying phase 
dynamics. 



VII. CONCLUSIONS 

We have discussed various phenomena in driven- 
dissipative many-body boson systems governed by many- 
body master equations with a dissipative zero mode. 
This includes both a semi-quantitative characterization 
of the steady state phase diagram resulting from the 
competition of the driven-dissipative dynamics with a 
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Hamiltonian, as well as a first characterization of meso- 
scopic nonlinear dissipative equilibration dynamics fol- 
lowing phase quenches in the Liouville operator. For this 
purpose, we have developed a generalized Gutzwiller ap- 
proach which accounts for density matrices correspond- 
ing to mixed states. We have introduced a low density 
limit which exhibits all the qualitative features of the 
phase diagram and is amenable to an analytic discussion 
due to the decoupling of a finite subset of correlation 
functions. 

Several features are direct consequences of the fact 
that our system is far out of thermodynamic equilibrium. 
Most prominently, the phase transition in steady state 
shares features of a quantum phase transition in that it 
is interaction driven, and of a classical phase transition 
in that the ordered phase terminates in a thermal state, 
in stark contrast to dissipative quantum phase transi- 
tions in global thermodynamic equilibrium, where the 
constraint of zero temperature inhibits such behavior. 
Furthermore, in the limit of vanishing dissipation, the 
nonequilibrium phase diagram for the steady state does 
not connect smoothly to the phase diagram of the equi- 
librium Bose-Hubbard model - this point is a manifesta- 
tion of the non-commutativity of the limits k — >• and 
t oo. Another unconventional signature is the exis- 
tence of a dynamical instability, which manifests itself in 
a spontaneous breaking of translational symmetry. Fi- 
nally, we have identified interesting nonlinear dissipative 
dynamics following a quench from one dark state to the 
other in mesoscopic systems, which hints at an interest- 
ing interplay of fluctuation and dissipation for effective 
macroscopic topological variables [27l |40] . 

We believe that the driven-dissipative many-body sys- 
tems bear substantial potential for intriguing many-body 
physics far away from thermodynamic equilibrium. In 
particular, interesting questions concern the exploration 
of the critical behavior close to the phase transition and 
a more analytical understanding of the mesoscopic dissi- 
pative dynamics in terms of collective variables. 
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Appendix A: Derivation of the mean field EOM 

In this section we compute the right-hand side of 
Eq. (J8|. Let us start from the local interaction term 
in the Hamiltonian 

v 



= Tr ^E Pt"[Unv(n v - I), p & >] . (Al) 
v v±v 

In the summation above we separate the term with £' = £ 
from the others. In the isolated term the trace applies 
only to the tensor product of density matrices and the 
commutator remains unaltered. In the remaining terms, 
instead, the trace applies also to the commutator part 
and hence the contribution vanishes. Hence, we obtain 

TrM^2 u Mnt' ~ 1)>P] = [Uniint - l),p e ] . (A2) 

v 

The contribution of the chemical potential, similarly, is 
the commutator with the term —pfi£. To treat a generic 
nonlocal term Oa> in the Hamiltonian or in the Lind- 
blad structure we introduce the decomposition Ow = 
E rs l r d> A \ ® Bf,. The hopping term b\b t + b\,b t on the 
^the site reads in this form 

= <*ra, A\ = (bi, b\) r , B\ = (b\M)r , (A3) 

with £ and £' nearest neighbors and 5 rs the Kronecker 
symbol. The trace of the commutator gives 

V rs 

= Tt^E Pt"'E2'yt>e» A t'® B e»>Pi'®Pt"] 

V £">jte',e" rr' 

(A4) 

Breaking the sum over £' as described above we obtain 

tt^[-j^SjW] = - J E + Pt] > 

IV (V\£) 

(A5) 

where (£'\£) indicates the sum over all the nearest neigh- 
bors £' of £. Finally, the contribution of the Hamiltonian 
is Tr ^t[H, p(t)] = [hi(i), pe(i)], with the effective locally- 
acting Hamiltonian 

hi = -J ({be>)b\ + {b\ f )be) + ^Uh £ (he - 1) - phi . 

(V\£) 

(A6) 

In the Liouvillian term J2(ij) ^^tfaijpclj ~ c lj c ijP ~ 
pc\jCij] in Eq. |8| we insert the identity 

1 = l — Sj£—5i£ + 5j£+5i£ + 5i£5j£ = (l — 5i£)(l — 5j£)+5i£+5j£ , 

(A7) 

that holds since j ^ i and so SjiSji — for every value 
of £. The summation breaks into three parts 

f E +E^ + E^) Tr ^[---]- (as) 

\(ij),i^£,i^£ (i\£) (j\£) j 

The first term gives zero, since the expression Tr^[. . .] in 
this case is equivalent to p^Tr[. . .] that vanishes for every 
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couple of i and j because of the form of the Lindblad 
term. Using that cn> = — Q/^, the two remaining terms 
give 



2 J2 ^M 2c i' 

(£>\£) 



£ — C £/£ C£'£P — pC £ , £ C£'£ 



(A9) 



To proceed further, we assume a cubic square lattice ge- 
ometry, with the coefficients 4>u> that vanish on each 
plane orthogonal to the axial direction x and are equal to 
a constant <fi along x. Then the phase difference between 
adjacent sites in the lattice is 4>w — vfy, with a = £' — t 
along the x axis and a = in the orthogonal plane. The 
jump operators in the decomposition above read 

lw = diag(-l, -e^, e-^t 1) , (A10) 



A\ = (l,6j,6£,n € ) r , B\ = (n € ,6 € ,6j,l) r . (All) 



The contribution of a single neighbor in Eq. (A9) is 
YpA\ Pl Af - A*A\ Pt - Pi AfA\)Y 



(A12) 



with the matrix of correlation functions 



/ (n 2 t )i 



(n e b e ) e e i °<f> 



(bine) i 



-(n e bl) e e 
~(ne)e 



(ne) 
-(bf)ee 

~(b\ 



-ilad 



-(b l n l ) l e i °* -(ne)e \ 
-(bj) f e i2 ^ 

(ne). 
(be)ie iu > 



1 



A 



~(be)t 
(b\)ee 
1 

(A13) 

The contribution of a — 1 and a = — 1 appears once, 
while cr = appears z — 2 times, where z is the coordi- 
nation number of the lattice (4 for the two-dimensional 
square, 6 for the three-dimensional cube). 



Appendix B: Linearization of the EOMs around the 
thermal state 



We begin the analysis with the definition of the den- 
sity matrix pi(t) = pf^ + 5pi(t) linearized around the 
thermal state and we keep only the first order in the en- 
tries of the perturbation 5 pi. The density matrix enters 
both the Hamiltonian hi = hff^ + 5hi and the Liouvillian 

Ci = df^ + SCi in the computation of the correlations. 
To the first order, the EOM for the perturbation matrix 
reads 



the Liouvillian. The terms in Eq. (Bl) read 



-bib 



Sh e (t) 



(B2) 



JbeJ2 s (bl +a )-JblJ2 S (^+-) 



pA^A\\ 



-,(0) Tt 



=±1 



8C e [p] = «5^[2AJpAf -A s ?A\p-pA s ?A\\ £ 5TJ« 



= ±1 



where we have reserved the possibility for the perturba- 
tion to be inhomogeneous in real space. The unperturbed 
matrix of the zeroth-order correlations, computed on the 
thermal state, reads 



/2n( 



V 



n + 1/2) 



-n 






71+ 1 





(B3) 



The matrix of the first-order correlations is 



STt. 



\ 



5(nf) 
5(nebe) 
-5(b e h e )* 
-5(he) 



5(hebe)* 
5 (he) 

S(b 2 e )* 
S(be)* 



-8 (b e ne) 

-Hi> 2 e) 

6(he) 
5(be) 



-S(n e ) \ 
-5(b e ) 

5(kr 





(B4) 

with the averages S(. . .) = Tr[. . . 5 pi]. We verified with 
the numerical integrations of the nonlinear EOM ([sj) that 
5 (of) is much smaller than the other correlation functions 
and hence can be neglected. 

The structure of the EO M for 5 pi is determined by the 
right-hand side of Eq. (Bl ), that is tridiagonal in the Fock 
space. It follows that only the entries of 5 pi on the upper, 
central, and lower diagonal need to be taken different 
from zero. Moreover, the equations for the entries in each 
diagonal constitute a closed system, so that we have only 
to consider the variables x v = bpi-,v,v-i-> where v > is 
the index in the Fock space. 

We derive a set of rules to project the EOM onto the 
lower diagonal. For tridiagonal matrices (such as 5 pi) it 
holds that 



P 

pS+S 

pbWb 
bpti 



— )► 
— )► 



Xu 

(y - l)x v 
{y-\fx v 



bWbp 
b%tib 
Ppb 



vx v 
v 2 x v 



and for diagonal matrices (such as 



— )> 



v(y — 1)q 



(B5) 

f h ') it also holds that 



dt6pt(t)rt[h?\5 Pe {t)]-C? ) [6pt(t)] = -i[8he,p e th) ]+8Ce[p 

(Bl) 

where the left-hand side is linear in the perturbation ma- 
trix, while the right-hand side acts as an effective driv- 
ing. The right-hand side depends on the perturbation 
density matrices Spi±i on the neighboring sites of £ via 
the correlation functions that enter the Hamiltonian and 



*(0)r 



(th) 



tiptib 
tibpti 

pbW 

pftbb 



— > 

— > 
— > 






y/v(y — l)x 
\[vvx v 

y/vVX v 





v-l 



&v 

bppb 
tibpb 

frfrbp 

bp bp 



yJVX v 




yjv(y - l)x v - 
. 



(B6) 
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We apply these rules to the EOM (Bl) for the pertur- 



bation matrix and we obtain the EOM for the variables 



2k 



d t x v = - [-i/(4 + 8n + iU/hz) + iU / n)x v 

+ 2ny/v(v - V)x v -\ + 2(n + \)^v(y + 1)^+1 

+ v^(<^-i> + <^ + i»(/4 h -i,-i + ft) 

+ ^[-«y(S/_i) (2i/ -»j/k) 



5<S m )(2i/ -*J/k) 

2((5(S^_in^_i) + 5(^ + in^+i))] 
Jth) X 



, (th) Jth) 



(B7) 



We see that in the EOMs the subdiagonal entries of 
the perturbation matrix on the site i are coupled to some 
correlation functions computed on the neighboring sites. 
To allow for any possible spatial configuration of the cor- 
relations makes the problem analytically unmanageable. 
Moreover, to compute the linearized equations on a fi- 
nite number of lattice sites induces a discretization of 
the Brillouin zone such that the numerical effort to solve 
the linear equations scales quadratically with the size of 
the lattice. At this point it is necessary to introduce 
some physical hypothesis on the form of the instability, 
which allows us to consider only a finite number of de- 
grees of freedom and to easily scan the Brillouin zone of 
the lattice. Here we make the hypothesis 



S(b e ) 
6(b e ,n e >) 



Sibtie-W-Q , 
5 (b e ne) e-^'-^ . 



(B8) 



Usin g E q. ( |B8[ ) we reduce the set of variables in the 
EOM ([B7|) to the sole {x u }. The EOM is rewritten in 
the form d t x v = J2 U , M u yx u >, where 



2n\Jv{y- 1)^/^-1 + i[-z/(4 + 8n + 



K 



.u 



+ i—]S v y + 2(n + l)yjv{v + 1)<^>+i 



.J 



+ v z/z/[cos (0o)2(i/ / — v) + i— cos ^oj 

K 

x (o {0) -o {0) ) 

X \Pi-v-\,v-\ Pt-v,v) 

+ 2cos(<^vW(^^ : 
as can be verified with the formal substitutions 

x„ -> b v y, 5(bi) ->> \[\ 



(B9) 



V, S(b £ h £ ) z/Vz/' . (BIO) 

The general solution of the EOM is given by x v ~ e A *, 
with A an eigenvalue of the matrix M. If Re [A] < for 
any A, the perturbation 5 pi damps out and the thermal 
state is stable. On the contrary, if there is at least one 
eigenvalue A with a positive real part 7 = Re [A] the ther- 
mal state is unstable. The corresponding eigenvector is 
the unstable mode in the channel labelled by 0q. To 



diagonalize the matrix M it is necessary to introduce a 
truncation z/ max in the set of variables {x u }, correspond- 
ing to a truncation in the Fock space of the system. We 
use z/max = 100 and we have verified in some typical cases 
that the higher part of the spectrum remains unchanged 
when z/ max is increased. 



Appendix C: Linearization of the EOMs around a 
generic mean field steady state 

In this section we outline the computation of the sta- 
bility spectrum of a generic mean field homogeneous 
steady state with local density matrix p^. This method 



is used in Sec. Ill C to compute the phase border be- 
tween the stable and the unstable condensate region. 
The mean field EOMs (p| can be rewritten in the form 
dtpi(t) = Q[pi(t), pi±i(t)] : with Q a bilinear operator. 
Let us characterize now the action of the operator Q on 
the steady-state p^(too) of the system. We require the 
steady state to be such that d t \p£ ;u y(t oo)\ = 0, where v 
and v' are indices in the Fock space of the t\h site. It 
is easy to show that a sufficient and necessary condition 
for this to happen is d t p£ ;u y(too) = —iQw'Pi-^y, with 
£l vv > G R. We also require that the absolute value of a 
generic correlation function (b\ p b^) is constant in time, 
and hence d t (bl P b^) = —iA pq (bl P b^). The stationarity of 
the absolute values of the correlation functions can be 
easily checked on the numerical solution of the equations 
of motion. 

Still, we cannot forbid a residual phase precession both 
in the off-diagonal elements of the density matrix and in 
the correlation functions. The first point raises a problem 
for the linearization of the equations of motion because 
we do not obtain numerically a state that nullifies exactly 
the bilinear operator Q, i.e. we do not have a good start- 
ing point for the instability analysis. It turns out that a 
chemical potential p can always be fixed into the operator 
Q (that we denote after the choice is made) such that 
Q/j,[P£(toc), Pe± i(too )] = exactly. To explain why this is 
the case, we consider the explicit form of the correlation 
functions <S| P S|> = ^ f(q, v + q)f(p, v)pt- v + q ,v+ v , with 
f(Pi v) = (v -\- p\b^ p \v) . From the conditions above it fol- 
lows that £}„+ p ^+ q = A pq and hence ft can only depend 
on the difference of its arguments and be of the form 
ftp- q . Now we may argue that Q should have the linear 
form Q x (p — q) since we want the phase frequency of 
(of) be twice that of (bi) in the condensate limit. Finally, 
Q should be a linear function of time if we require that 
we can remove the phase precession by a frame transfor- 
mation of the basis vectors. The result is that, in the 
steady state, the residual time-evolution of the density 
matrix is of the form dtpe-^y = —iuj{y — v')pt,vy. This 
is precisely the contribution to the equation of motion 
produced by a term uoh in the Hamiltonian and can thus 
be eliminated by the appropriate choice of a chemical 
potential p = — uo. 
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With such choice of chemical potential, the equation 
of motion for the fluctuation 5pi(t) — pi(t) — poo of 
the density matrix reads d t Spi(t) = Q^Jpoo? 5pe±i(t)] + 



Appendix D: EOM of the correlation functions 



Qn[dpi(t), Poo}- Then, as we point out in Sec. IIIC 



well, some hypotheses on the spatial dependenceofthe In this section we provide the EOMs of the time- 



unstable mode has to be made to relate the density ma- 
trices on different sites as 5p£±\ = T^[5p£ : (j)o] : where 
we restrict to a one-dimensional geometry. According 
to the results of the numerical solution of the EOM ([8]), 
here we consider the stability of the modes defined by 
[T^ [Spi, 0o]W' = cos (0o)[/^W'j that correspond to 
the generation of a standing wave in the lattice with 
wavelength 27r/0o. The EOM for the fluctuation ma- 
trix then reads Spi(t) = L[5p£, <po] with the linear opera- 
tor L given by Q>^ )r] o (T(+>[., fa] + T<->[., O ]) + 
pi(to)]. The existence of a positive real eigenvalue 
of £[-,</>o] means that the steady state p^ is unstable in 
the particular mode labeled by 



dependent correlation functions c^qj = (b\ p bf) that fol- 



low from the mean field EOM (|8j) for the density ma- 
trix. For definiteness, in the following we assume a one- 
dimensional geometry so that the neighbors {£'\£) of the 
sites £ have indices £ + 1 and £ — 1. The EOM of a generic 



correlation function is dt{b p b\) — ^[b P b^dtpi{t)]. Sub- 
stituting the time derivative of the density matrix p£ one 
obtains products of averages of operators on the site £ by 
the mean fields ( 11 ) on the neighboring sites. Normal or- 



dering the operators acting on the site £, one obtains the 
following expression, which is linear in the correlation 
functions on the site £ but contains products between 
correlation functions on neighboring sites 



dtc p , q ,e = + [i(U/2)(p(p - 1) - q(q - 1)) - 2«((p - q) 2 +p + q) - ip(p - q)] c p , q/ 
+iU (p - q)c p+li q +li £ + 2Kpq(c lyl: £- 1 + Ci,i,£+i)c p -i, g -i,£ 

+ [-iJp(c lt0 ,£-i + ci >0 ,£+i) + 2kp(1 - q)(e i(f) c ly0 ,e-i + e" <0 ci jO ,£+i) + 2«p(e i0 c 2 ,i,£_i + e _i</) c 2 ,i,£+i)] Cp-i iqi i 
+ [iJqicoxi-i + c o,M+i) + 2Kq(l - p)(e"^c ,i^-i + e z0 c o ,i,£ + i) + 2«g(e~^ci >2 ,£_i + e^ci >2 ,£+i)] Cp,g-M 
-2^(e*^ci j0 ,£-i + e~*^ci j o,£+i)c p> g + i > £ - 2pn{e~^ c^ X t-\ + e^c ,i^ + i)c p+ i 5(? ^ 

- l)(e 2 ^c 2 , ,^-i + e" 2i</) c 2 ,o,£+i)c p _ 2 , (? ,£ + - l)(e _2 ^co >2j £_i + e^co^+i)^-^ , (Dl) 



where it is understood that c p ^£ — whenever p,q < 
and co,o,£ = 1 by the unitarity of the trace of p£. The 
above expression shows that the EOMs of the correla- 
tion functions are organized as an infinite hierarchy, since 
c Py q y £ is also coupled to c p+ i jP+ i^ (due to the local inter- 
action) and to Cp^+i^, Cp+i^qj (due to the dissipation). 
Hence we do not find a subset of EOMs that decouples 
exactly from the others and that can be solved in a finite 
number of steps. 

In Sec. (ITvT) we use the different set of "connected" 



correlation functions (56j p 56|), built with the fluctuation 
operator 5b £ = b£—ip£. We remark that we always work in 
the Schrodinger picture here: the operators 5b t are con- 
stant in time and depend parametrically on the choice 
of the quantity ip£. This set is totally equivalent from 
a mathematical point of view, but is motivated from a 
physical perspective if we interpret ip£ as the order pa- 
rameter of a Bose-Einstein condensate and (11) holds. 



{5b £) = 0. We reiterate that the vanishing of the first- 
order correlation takes place only when the average is 
taken with respect to some density matrix pt(to) such 
that (11) is valid. Once we let the density matrix evolve 
in time, as we do in Sec. [v[ the operators 5b £ do not 
change but the average Tr[5b£p£(t)] = Tr[b£p£(t)] — ip£ is 
in general different from zero. 

To obtain the EOM for the correlation functions 5b £ 
we expand the product of \)£ operators and ip£ c-numbers 
in the average, obtaining a sum of correlation functions 
c P: q y £. Then we substitute the EOM ( |D1| ) of each c VA ^ 
and again we expand the products of operators in the 
averages in terms of 5b £ and ip£. This approach does not 
yield the general EOM, analogous to (Dl), because each 



The new set of correlation functions hence describes the 
fluctuations about the condensate and, by construction 



choice of p, q produces a different set of correlations, but 
can be applied straightforwardly. The EOMs for the con- 
nected correlation functions are organized as an infinite 
hierarchy as well, and cannot be reduced analytically to 
a closed subset. 
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